function ivp

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

% clear all previous variables and plots
clear *
clf

% time interval
tmax=20;

%  initial values
y10=0.5; y20=0; 
y0=[y10 y20];

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);  titleX='ode45';
%[t,y] = ode23s(@rhs,[0 tmax],y0);  titleX='ode23s';

%
%  2-D plot:  plot solutions as function of t
%
hold on
plot(t,y(:,1),'r')
plot(t,y(:,2),'--b')

% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')

% command to put legend into plot
legend(' y_1',' y_2',4);

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 

hold off

%
%  3-D plot: plot solution vector
%
figure
hold on
plot3(t,y(:,1),y(:,2),'b')
view(50,15)
grid on
box on

% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('y_1-axis','FontSize',14,'FontWeight','bold')
zlabel('y_2-axis','FontSize',14,'FontWeight','bold')

% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
hold off


%  define f(t,y)
function dy=rhs(t,y)
dy=zeros(2,1);
a=8e-4; ep=4e-2; f=1/4;
dy(1) = ( y(1) - y(1)*y(1) + f*(a-y(1))*y(2)/(a+y(1)) )/ep;
dy(2) = y(1) - y(2);













